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SECTION  I 


INTRODUCTION 


Researchers  in  aerodynamic  analysis  and  design  have  turned  to  numerical 
techniques,  with  the  advent  of  the  supercomputer  (Cyber  205,  Cray  X-MP,  Cray 
2,  etc.),  to  solve  for  the  fluid  flow  about  weapon/store  configurations. 
Current  aerodynamic  research  at  the  Air  Force  Armament  Laboratory  (AFATL)  is 
aimed  at  solving  the  equations  that  govern  fluid  flow  problems  using  a 
variety  of  approximation  techniques.  The  governing  partial  differential 
equations  (PDE)  form  a  nonlinear  system  which  must  be  solved  for  the  unknown 
pressures,  densities,  temperatures,  and  velocities  to  yield  the  aerodynamic 
characteristics  for  a  given  weapon/store  conf iguration  at  specific  flight 
conditions . 

To  obtain  a  thorough  approximation  for  the  flow  field  about  a  configu¬ 
ration,  we  must  solve  the  complete  Navier-Stokes  equations.  However,  due  to 
limitations  placed  on  researchers  by  current  computer  systems  (time  and 
storage),  certain  simplifying  assumptions  must  be  made  to  obtain  results. 

By  assuming  an  inviscid,  adiabatic  flow  field  (dropping  viscous  and  heat 
transfer  terms),  we  obtain  the  Euler  equations.  Results  obtained  from  a 
solution  of  the  Euler  equations  are  particularly  useful  in  preliminary 
design  work  where  information  on  pressure  alone  is  desired.  These  equations 
are  also  of  interest  because  they  incorporate  many  major  fluid  dynamics  ele¬ 
ments  such  as  internal  discontinuities  (shock  waves  and  contact  surfaces). 

The  Euler  equations  govern  the  motion  of  inviscid,  non -heated  gas  and  have 
different  numerical  characteristics  in  different  flow  regimes.  For  steady 
problems,  the  equations  are  elliptic  in  subsonic  flow  and  hyperbolic  in 
supersonic  flow  (Reference  1  ). 

Research  at  AFATL  is  aimed  at  solving  the  three  dimensional  Euler 
equations  to  approximate  the  flow  about  arbitrarily-shaped  weapon/store 
configurations.  This  research  is  accomplished  using  various  implicit  and 
explicit  finite-difference  and  finite-volume  techniques.  Currently  this 
work  has  lead  to  an  analysis  of  two  explicit,  finite-volume  procedures  for 
solving  the  Euler  equations.  To  help  in  this  analysis  of  the  numerical 
characteristics  inherent  to  the  two  approaches,  a  single  equation  serving 
as  a  numerical  analog  to  the  Euler  equations  has  been  found.  The  inviscid 
Burgers  equation  (nonlinear  wave  equation)  serves  as  this  simple  nonlinear 
analog  to  aid  in  our  understanding  of  these  techniques  (Appendix  A) 

(Reference  1). 


SECTION  II 


UPWIND  SCHEME 


Several  numerical  techniques  have  been  developed  that  will  solve  the 
partial  differential  equations  that  govern  fluid  flow,  heat  transfer,  and 
combustion  problems.  These  methods  can  be  either  explicit  or  implicit, 
central  difference  or  upwind,  single  or  multi-step,  and  first-  or  second- 
order  accurate. 

An  explicit,  second-order,  one-sided,  or  upwind  difference  scheme  for 
the  numerical  solution  of  hyperbolic  systems  has  been  developed  by  Warming 
and  Beam  (Reference  2).  There  are  several  advantages  to  the  use  of  the 
Upwind  schemes. 

(1)  One-sided  schemes  are  often  desirable  along  both  fixed  external  ' 
boundaries  and  along  moving  internal  boundaries  (such  as  shocks),  where  a 
spatially  centered  scheme  would  require  one  or  more  points  inside  or  across 
the  boundary. 

(2)  An  explicit,  second-order  upwind  scheme  can  have  twice  the 
stability  bound  of  a  symmetric  scheme  using  the  same  number  of  spatial  grid 
points . 


(3)  The  dissipative-dispersive  properties  of  an  upwind  scheme  are 
superior  to  those  of  a  symmetric  scheme. 


(4)  Switching  schemes  across  a  discontinuity  can  reduce  the  spurious 
oscillations  usually  associated  with  a  second -order  accurate  shock-capturing 
technique.  Warming  and  Beam's  goal  was  to  develop  a  hybrid  scheme  which 
exploits  the  advantages  of  a  second-order  upwind  scheme  for  aerodynamic 
flows.  This  multi-step  procedure  applied  to  the  inviscid  Burgers  equation 
is  shown  as  follows: 
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A  thorough  error  analysis  has  been  performed  on  the  Upwind  scheme  to 
aid  in  our  understanding  of  the  numerical  characteristics  inherent  to  this 
technique  (Appendix  B).  This  analysis  includes  consistency,  convergence, 
numerical  stability,  phase  and  dispersion  error,  and  artificial  dissipa¬ 
tion  (Reference  3). 


SECTION  III 


FINITE-VOLUME  APPROACHES 


For  our  research  in  computational  aerodynamics,  we  have  settled  on  the 
finite-volume  (FV)  approach,  as  opposed  to  the  finite-difference  (FD) 
method,  to  solve  for  the  physics  of  the  flow  about  a  weapon/store  configura¬ 
tion.  Fundamentally,  the  primary  difference  between  the  two  methods  is  that 
for  the  FV  approach  we  solve  for  the  flux  at  the  face  of  a  cell;  the  FD 
approach  solves  for  the  flux  at  the  center  of  the  cell  (Figure  1 ).  To 
accomplish  this  FV  technique  we  must  extrapolate  either  the  flux  or  the 
dependent  variable  from  the  center  of  the  cell  to  the  face  of  the  cell.  The 
advantage  to  using  the  FV  method  is  that,  inherent  to  the  approach,  the 
conservative  property  of  the  PDE  is  fully  maintained. 

1 .  EXTRAPOLATION  TECHNIQUES 

As  mentioned  above,  there  are  two  extrapolation  techniques  that  must  be 
studied  to  determine  which  yields  the  optimum  results.  The  flux  term  is 
U(i)2/2  where  the  dependent  variable  is  U(i).  For  the  upwind  scheme, 
extrapolating  the  flux  yields 

2  (U(i)2/1)  -  (U(i-1)  2/2 )  or  (3) 

2  ( U(i+1 )  2/2 )  -  ( U(i+2)  2/2 ) 

depending  on  the  direction  of  the  flow  of  information.  When  extrapolating 
the  dependent  variable  for  the  same  inviscid  Burgers  equation,  we  get; 

( 2U(i )  -  U(i-1))  2/2  oi  (4) 

(2U(i+1 )  -  U(i+2) )  2/2 , 

again,  depending  on  the  direction  of  the  flow  information.  In  this  way,  the 
dependent  variable  extrapolation  technique  yields 

2U(i ) 2  -2U(i)U(i-1)  +  U(i-1 ) 2/2  or  (5) 

2U(i+1 ) 2  -2U(i+1 )U(i+2)  +  U(i+2)2/2  . 

This  result  shows  an  extra  term,  -  2U(i+1 )U(i+2)  ,  (as  compared  to  Equation 
3)  which  will  have  some  effects  that  are  inherent  to  this  type  of  numerical 
method . 

2.  DIFFERENCING  TECHNIQUES 

The  flow  solvers  (Euler  codes),  being  examined  by  our  research,  use  two 
types  of  differencing  techniques  to  solve  for  the  flux  at  the  face  of  a 
cell;  both  of  these  techniques  have  been  developed  by  Whitfield  (References 
4,  5,  and  6). 


The  first  differencing  technique  (Figure  2)  employs  a  dependent  variable 
averaging  (DVA)  approach  to  determine  the  direction  of  the  flow  of  informa¬ 
tion  across  a  cell  face.  For  a  specific  cell  face,  the  value  of  the  depen¬ 
dent  variable  (U)  is  averaged  from  both  sides  of  the  cell.  This  yields  a 
value  for  the  dependent  variable  and,  depending  on  the  sign  of  U,  is  used  to 
extrapolate  (using  one  of  the  techniques  discussed  above)  from  one  side  of 
the  cell  face  or  the  other,  to  obtain  a  value  for  the  flux  at  the  face  of 
the  cell. 

The  second  differencing  technique  (Figure  3)  employs  a  dual  dependent 
variable  technique  ( DDV)  in  which  the  value  for  the  dependent  variable  (U) 
is  determined  for  both  sides  of  the  cell  face  and,  depending  on  the  sign  of 
U,  can  utilize  both  values  of  the  dependent  variable  if  the  direction  of  the 
flow  of  information  is  toward  the  cell  face  from  both  directions.  The 
direction  of  the  flow  of  information  determines  whether  the  value  of  the 
dependent  variable,  CJ,  is  used  from  one  side  of  the  cell  face  (or  the  other) 
or  from  both  sides,  to  extrapolate  the  flux  to  the  cell  face. 


SECTION  IV 


NUMERICAL  RESULTS 


For  this  analysis,  both  the  extrapolation  and  differencing  techniques  are 
examined  using  the  inviscid  Burgers  equation  to  model  the  propagation  of  a 
wave  in  time.  The  numerical  analysis  used  here  to  study  the  inviscid 
Burgers  equation  should  yield  second-order  accurate  results.  This  level  of 
accuracy  is  typified  by  dispersion  or  ringing  effects  which  overshoot  the 
actual,  physical  results. 

The  first  phase  of  our  study  looks  at  forcing  a  wave  to  propagate  in  only 
one  direction.  This  allows  us  to  better  examine  the  differences  between  the 
two  extrapolation  techniques;  since  the  wave  is  moving  in  only  one  direc¬ 
tion,  the  differencing  techniques  are  essentially  the  same.  Figure  4a  shows 
the  effects  of  using  the  DVA  technique  with  dependent  variable  extrapola¬ 
tion.  The  results  are  atypical,  for  a  second-order  accurate  solution 
method,  in  that  no  dispersive  effects  (ringing)  are  apparent  downwind  of  the 
wave.  Figure  4b  gives  the  results  for  the  DVA  technique  using  flux  extra¬ 
polation.  These  results  show  the  characteristic  dispersive  effects  yielded 
by  a  second-order  scheme  in  which  the  numerical  solution  overshoots  the 
exact  solution  on  the  downwind  side  of  the  shock.  Figure  4c  shows  results 
for  the  DDV  technique  with  dependent  variable  extrapolation  in  which  we 
again  observe  what  appears  to  be  dissipative  characteristics  to  a  second- 
order  scheme.  In  Figure  4d  the  results  are  given  for  the  DDV  technique 
using  flux  extrapolation.  Once  again  the  flux  extrapolation  approach  yields 
typical  second-order  results  with  ringing  effects.  To  better  study  the 
effects  of  the  extrapolation  approaches,  a  Courant  number  of  1.50  was  used 
since  at  a  Courant  number  of  either  1.00  or  2.00  the  numerical  technique 
yields  the  exact  solution  (by  satisfying  the  shift  condition  in  Equation 
B.9  )  (Reference  1  ). 

The  second  phase  of  our  examination  studies  the  effects  of  the  two  differ¬ 
encing  approaches  (DDV  and  DVA)  by  forcing  two  waves  to  meet  at  a  cell  face 
of  equal  velocities  (magnitudes).  A  Courant  number  of  2.00  is  employed  so 
that  the  effects  of  the  two  extrapolation  techniques  are  negated  (by  satis¬ 
fying  the  shift  condition  in  Equation  B.9)  (Reference).  Figure  5a  shows  the 
results  for  the  DVA  technique  using  dependent  variable  extrapolation  which 
yields  the  exact  solution  to  the  mathematical  model.  Similar  results  were 
yielded  in  Figure  5b  for  the  DVA  technique  using  flux  extrapolation.  Figure 
5c  yields  the  results  for  the  DDV  technique  with  dependent  variable  extrapo¬ 
lation  and  shows  the  exact  solution,  as  does  Figure  5d  for  the  DDV  technique 
using  flux  extrapolation.  For  this  simple  model  both  differencing  tech¬ 
niques  (DDV  and  DVA)  yield  the  same  results;  therefore,  a  more  complicated 
test  must  be  accomplished  to  better  understand  the  limitations  of  the 
approaches . 

For  the  final  phase  of  this  investigation,  all  aspects  of  the  problem  are 
examined  by  forcing  two  waves  to  approach  each  other  at  unequal  velocities 
(magnitudes).  Due  to  the  unequal  velocities,  U  ,  the  effective  Courant 
number  changes  for  each  direction.  Figure  6a  shows  the  effects  for  the 
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Figure  4c.  Burgers  Equation  Solution  -  \>  =  1 . 50 
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igure  5a.  Burgers  Equation  Solution  -  u ~2 . 00 
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DDV  technique  with  dependent  variable  extrapolation.  The  results  show  the 
atypical  (dissipative)  solution  yielded  by  the  dependent  variable  extrapola¬ 
tion.  As  the  waves  meet,  the  wave  with  the  greater  magnitude  runs  over  the 
lesser  wave  at  an  average  wave  speed.  In  Figure  6b  we  obtain  the  results 
for  the  DVA  technique  using  flux  extrapolation.  This  approach  yields  the 
expected  second-order  results  and  also  shows  the  greater  magnitude  wave 
running  over  the  lesser  wave  at  an  average  (or  deduced)  velocity.  Figure  6c 
shows  results,  similar  to  those  in  Figure  6a,  for  the  DDV  techniques  using 
the  dependent  variable  extrapolation  approach.  As  in  Figure  6a,  we  observe 
the  dissipative  effects  on  the  solution  and  the  propagation  of  the  larger 
wave  at  an  average  speed  after  the  collision.  Figure  6d  shows  the  results 
for  the  DDV  technique  with  flux  extrapolation.  As  in  Figure  6b,  the 
expected  dispersive  effects  are  observed  and,  again,  the  greater  wave  moves 
at  an  average  velocity  after  the  collision. 

A  further  analysis  was  performed  to  attempt  to  find  weaknesses  in  the  dif¬ 
ferent  approaches.  In  this  case  the  effective  Courant  number  was  doubled 
thereby  pushing  the  CFL  condition  towards  its  limit  of  2.00.  Figure  6e 
shows  results,  similar  to  those  in  Figure  6a,  for  the  DVA  technique  with 
dependent  variable  extrapolation;  however,  the  propagation  of  the  wave  after 
the  collision  has  twice  the  velocity.  In  Figure  6f  the  results  are,  again, 
similar  to  those  in  Figure  6b  with  the  exception  of  the  wave  propogation 
being  at  twice  the  speed.  Figure  6g  shows  the  only  major  discrepancy  in  our 
analysis.  The  DDV  technique  with  dependent  variable  extrapolation  shows 
good  results  until  the  meeting  point  for  the  unequal  waves.  At  the  colli¬ 
sion  due  to  diffusive  effects,  the  effective  Courant  number  goes  well  above 
the  stability  limit  of  2.00  and  the  solution  diverges.  However,  when  the 
DDV  technique  is  applied  with  flux  extrapolation  (Figure  6h),  we  again 
obtain  typical  results  similar  to  those  in  Figure  6d. 
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Figure  6b.  Burgers  Equation  Solution 
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Figure  6e.  Burgers  Equation  Solution 
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Figure  6h.  Burgers  Equation  Solution 


SFCTION  V 


CONCLUSION 


This  investigation  has  yielded  several  important  conclusions  dealing  with 
both  the  extrapolation  techniques  and  the  differencing  approaches. 

(1)  The  dependent  variable  extrapolation  technique,  used  for  both  the 
dual  dependent  variable  (DDV)  and  the  dependent  variable  averaging  (DVA) 
differencing  methods,  tends  to  negate  the  typical  dispersive  effects  found 
in  second-order  accurate  schemes.  This  is  primarily  due  to  the  extra  term 
found  when  applying  the  dependent  variable  extrapolation  to  the  inviscid 
Burgers  equation.  When  applying  the  flux  extrapolation  technique  to  both 
differencing  methods,  the  expected  dispersive  effects  are  obtained. 

(2)  No  significant  differences  were  apparent  when  comparing  the  differ¬ 
encing  techniques  (dual  dependent  variable  and  dependent  variable  averag¬ 
ing).  Both  techniques  closely  modeled  the  correct  mathematical  solution, 
when  using  the  flux  extrapolation  approach,  throughout  most  of  the  investi¬ 
gation.  The  only  inconsistency  in  the  study  occurred  when  the  DDV  differ¬ 
encing  approach  was  used  with  dependent  variable  extrapolation  (Figure  6g). 
In  general,  this  tends  to  show  that  the  DVA  approach  may  be  a  more  robust 
method  than  the  DDV  approach. 

Therefore,  the  recommended  approach  to  solving  the  inviscid  Burgers  equation 
for  this  examination  is  to  apply  the  dependent  variable  averaging  differ¬ 
encing  technique  with  flux  extrapolation.  This  method  will  also  be  directly 
applied  to  the  three-dimensional  Buler  equations  for  the  solution  of  the 
inviscid  flow  field  about  arbitrarily  shaped  weapon/store  configurations. 


APPENDIX  A 


INVISCID  BURGERS  EQUATION 

Burgers  equation  (1948)  can  serve  as  a  nonlinear  analog  of  the  fluid 
mechanics  equations.  This  single  equation  has  terms  that  closely  duplicate 
the  physical  properties  of  the  fluid  equations,  i.e.,  the  model  equation  has 
a  convective  term,  a  diffusive  or  dissipative  term,  and  a  time-dependent 
term.  (Reference  1 ) 

Ut  +  UUX  =  U^  (A.  1  ) 

Unsteady  term  Convective  term  Viscous  term 

Equation  (A.1)  is  parabolic  when  the  viscous  term  is  included  and  is  a  good 
model  for  the  boundary-layer  equations,  the  parabolized  Navier-Stokes  (PNS) 
equations  and  the  complete  Navier-Stokes  equations.  If  the  viscous  term  is 
neglected  (inviscid  Burgers  equation),  the  remaining  equation  is  composed  of 
the  unsteady  term  (Ut)  and  the  nonlinear  convection  term  (UUX).  The 
resulting  hyperbolic  equation 


°t  +  UUx  =  0;  (A. 2a) 

in  conservation  law  form 

Ut  +  ( U2/2)  =  0,  (A. 2b) 

may  be  considered  a  simple  analog  of  the  Euler  eqautions  for  the  flow  of  an 
inviscid  fluid.  The  analogy  can  be  drawn  due  to  the  fact  that  both  equa¬ 
tions  are  first-order,  hyperbolic,  quasi-linear,  partial  differential  equa¬ 
tions  and  both  model  discontinuities,  such  as  shocks,  in  the  flow  field 
(Reference  1). 

Burger  equation  is  a  partial  differential  equation  because  U=U(x,t)  and  Ufc 
and  Ux  are  both  used  in  the  equation.  It  is  a  first-order  equation  because 
only  first  derivatives  appear.  Equations  (A.2a)  and  (A.  2b)  are  nonlinear 
because  the  unknown  variable,  or  dependent  variable,  U,  is  multiplied  by 
itself  or  its  derivatives.  Burgers  equation  belongs  to  a  special  class  of 
nonlinear  equations  called  quasi-linear  equations.  A  quasi-linear  equation 
is  one  in  which  the  highest  order  derivatives  appear  to  the  first  power,  as 
is  most  easily  seen  in  Equation  (A.2b).  The  three-dimensional  Euler  equa¬ 
tions  of  inviscid  fluid  flow  are  another  example  of  quasi-linear  equations. 

To  provide  a  good  test  case  for  the  finite  volume  approaches,  we  need  to 
know  some  exact  solutions  to  Burgers  equation.  First,  solve  Equation  (A.2) 
with  the  general  initial  conditions, 

U(x, 0)  =  f (x) .  (A. 3) 

The  solution  to  the  problem  will  be  of  the  form  U=g(x,t).  Now  consider  a 
three-dimensional  space  with  x,  t,  and  U  coordinate  axes,  and  define 


The  surface  <P  (x,t,U)=0  in  the  three-dimensional  space  defines  the  solution 
to  the  problem.  In  other  words,  any  point  on  the  surface,  say  xQ,  tQ,  UQ  i 
a  point  of  the  solution  given  by  U{xQ,t0)=  UQ.  A  vector  perpendicular  to 
the  surface  is  given  by  grad  <J>  ■  (4>*x,  <J>y ,  4>z),  and  any  vector  perpendicular 
to  grad  <p  must  be  tangent  to  the  solution  surface.  If  (xQ,  tQ,  UQ)  is  a 
point  of  the  solution,  and,  for  infinitesimal  displacements,  (xQ+dx,  tQ+dt, 
UQ+du)  is  also  a  point  of  the  solution,  then 


grad  <j>  •  dx  =  0, 

or  <J>xdx+  <J>  tdt+  ‘frydu  =0.  , 

Using  Equation  (A. 4)  in  Equation  (A. 5) 

-gx-gt  +  du  =  0.  (A-6 

6.  Prom  Equation  (A. 2)  the  following  results  when  U=  g(x,t), 

gt  +  Ugx  =  °» 

from  which  we  can  tell  that  if  we  take  dx,  dt,  and  du  in  the  following 
ratios.  Equation  (A. 6)  will  be  satisfied: 

dt/1  =  dx/U  =  du/0.  (A. 1 

If  r  is  taken  as  a  parameter  along  the  solution  curve,  then  Equation  (A.7) 
can  be  written: 

dt/dr  a  1  ,  dx/dr  =  U  ,  du/dr  =  0;  (A. 8 

with  the  initial  conditions: 

t  =  0  ,  x  =  xQ  ,  U  =  f(xQ). 

Solving  these  equations  gives 

t(r)  =  r  ,  U(r)  =  constant  =  f(xQ)  and 


dx/dr  gives 


x(r)  =  rf(xQ)  +  constant  and  x(0)  =  constant  =  xQ. 


Therefore; 


x(r )  =  rf(xQ)  +  xQ 
x ( t )  =  tf(xQ)  +  xQ. 


For  an  example,  take 


1. )  For  xQ  >  1  (Figure  A— 1 ) ; 

x(t)  =  tf(xQ)  +  xQ  =  xQ 

U(r)  =  f(xQ)  =  0 

U(x,t)  =  0. 

2. )  For  xQ  <  0  (Figure  A-1); 

x(t)  =  t  +  xQ 

U(r)  -  f(xQ)  =  1 

U(t+xQ/t)  =  1 

since  xQ  =  x-t  when  U(x,t)  =  1. 

3. )  For  0  <  x  <  1  (Figure  A-1); 

xQ  =  ( x— t ) / ( 1  — t ) 

U(x, t)  =  ( 1 —x )/( 1— t) • 

This  method  of  characteristics  analysis  shows  the  Burgers  equation  capabil¬ 
ity  to  produce  shocks  or  discontinuities  in  the  flow  field.  This  capability 
allows  one  to  test  numerical  shock-capturing  methods  using  the  inviscid 
Burgers  equation  and  apply  them  to  the  three-dimensional  Euler  equations  for 
inviscid  fluid  flow. 


APPENDIX  B 


ANALYSIS  OF  NUMERICAL  TECHNIQUE 


This  section  of  the  analysis  examines  the  numerical  characteristics  of 
the  Warming-Beam  algorithm  for  the  general  finite-difference  approach 
(References  1  and  2).  Finite-difference  methods  involve  approximating  the 
continuous  domain  of  any  problem  by  a  discrete  domain  (grid)  and  approxi¬ 
mating  the  PDE's  governing  any  problem  by  one  or  more  algebraic  or  finite 
difference  equations  (FDE's).  The  total  error  in  the  solutions  of  FDE's  is 
made  up  of  discretization  error  and  stability  error.  Stability  error  is 
small  for  stable  FDE's  since  by  definition  disturbances  and  errors  cannot 
grow.  Therefore,  discretization  error  accounts  for  most  of  the  total  error. 
The  discretization  error  is  made  up  of  dissipation  and  dispersion  error. 

For  this  examination  the  governing  PDE  is  the  inviscid  Burgers 
equation: 


Ut  +  (U2/2)x  =  0  .  ( B. 1 ) 

To  determine  if  a  FDE  is  an  algebraic  analog  of  a  PDE, and  is  agood 
approximation  of  the  exact  solution  of  the  PDE,  we  must  involve  the  concepts 
of  consistency,  numerical  stability,  convergence,  phase  and  dispersion  error, 
and  artificial  dissipation  (Reference  3). 

1  .  LINEARIZATION 

For  the  purpose  of  linear  stability  theory,  we  employ  a  linearization 
method  to  the  nonlinear  PDE, 

Ut  +  UUX  =  0. 

The  resulting  linearized  PDE  is  the  convection  or  linear  wave  equation: 

Ut  +  CUX  =0,  ( B. 2 ) 

and,  for  continuity,  the  linear  equation  will  be  used  throughout  the  FDE 
analysis . 

Locally,  the  PDE  may  be  approximated  by  the  linear  PDE  with  constant 
coefficients  even  though  the  PDE  is  globally  nonlinear.  This  approach 
yields  reasonable  results;  however,  the  resulting  stability  criteria  (to  be 
discussed  in  more  detail  in  subsection  3)  is  necessary,  but  not  sufficient. 
The  resulting  Warming-Beam  upwind  scheme  for  the  linear  wave  equation  is  as 
follows : 


n+1  n  n  n 

Predictor:  U,  -U,  +  C  U,  -UL  t  =  0 


(B.3a) 


Corrector: 


Uj  -  1/2  (Uj  +  Oj  ) 
A  t 
2 


Ax 


(B.3b) 


2.  CONSISTENCY 


To  analyze  the  consistency  of  the  FDE  (Equation  B.3a,b)  used  to  model 
the  PDE  (Equation  B.2),  we  must  express  U^+  ,  U?+1,  in  terms  of 
Uj1  and  its  derivatives  by  using  a  Taylor  series  expansion: 


n  +  l  n  i  "At2  "  1  \t2 

U.  =  U.  +  U.  At  +  U.~ —  +  U.  — - —  + 

l  11  i2 !  l  3! 


2  3 

n  1  "  4At  8At 

Un  +  2A t  +  U.  +  U  .  + 

1  1  1  2 !  1  3 ! 


u"  -  a1  t  *  u"  -  u":  + 

i  l  i  2 !  1  3 ! 


(B.4a) 


(B.4b) 


(B.4c) 


The  reason  for  this  is  that  the  FDE  was  derived  by  applying  the  PDE  at  grid 
point  i  and  time  level  n.  Substituting  the  expansions  into  the  FDE  yields 


Ut  +  cux  =  0  (A  x2,  A  x  At,  At2), 


(B.5) 


therefore.  Equation  (B.5)  is  mathematically  equivalent  to  the  PDE  as  the 
grid  spacing  (a  x)  and  time-step  size  (A  t)  approach  zero. 

An  FDE  is  consistent  if  for  every  i  and  n; 

lim  FDE  =  PDE  (B.6) 

A  x  -►  0 
A  t  -+  0 

Consistency  measures  the  extent  to  which  an  FDE  approximates  a  PDE  in  some 
limiting  case.  A  consistency  analysis  was  performed  on  the  upwind  scheme 
and  it  was  found  that  the  upwind  scheme,  when  applied  to  the  convection 
equation,  is  unconditionally  consistent.  This  means  that  the 


n 


n 


(B.7) 


lim  (FDE)i  =  (PDE)i 
At  ->  0 
Ax  -*•  0 

n  n 

regardless  of  how  x  and  t  approach  zero.  Terms  such  as  A  tU^  or  a  xU^ 
can  be  added  to  the  FDE  and  the  resulting  FDE  will  remain  unconditionally 
consistent;  therefore,  there  are  an  infinite  number  of  unconditionally 
consistent  FDE's  for  a  given  PDE.  Consistency  alone,  however,  will  not 
guarantee  the  accuracy  of  the  solutions  of  the  FDE's  since  all  computations 
are  performed  by  using  finite  A  x  and  A  t. 

3.  STABILITY 

Numerical  stability  is  concerned  with  how  errors  propagate  as  the  solu¬ 
tion  is  advanced  in  a  time-like  variable,  and  is  a  concept  applicable  only 
to  parabolic  and  hyperbolic  PDE's.  An  FDE  is  stable  if  the  stability  error 
approaches  zero  or  does  not  grow.  A  given  FDE  may  yield  stable  or  unstable 
numerical  solutions  depending  upon  the  value  of  some  dimensionless  parameter 
(  v  =  c  A  t/  A  x).  The  need  to  obtain  stable  numerical  solutions  is 
critical  to  the  solution  of  a  given  PDE  since  only  stable  numerical  results 
have  a  chance  of  being  physically  meaningful.  Conditionally  stable  FDE's  are 
those  that  yield  stable  solutions  when  A  x  and  A  t  are  in  a  given  form  (  y  = 
cAt/Ax).  Unconditionally  stable  FDE's  are  those  that  give  stable  solutions 
for  any  A  x  and  A  t.  Unconditionally  unstable  FDE's  give  unstable  solutions 
for  every  Ax  and  At.  Typically,  stability  bounds  for  implicit  methods  are 
less  restrictive  than  those  for  explicit  methods  . 

There  are  many  different  mathematical  methods  for  analyzing  the 
numerical  stability  of  FDE's.  For  this  examination,  the  Von  Neumann  or 
Fourier  method  is  employed.  In  the  Fourier  method,  the  numerical  stability 
of  an  FDE  is  analyzed  by  introducing  a  disturbance  into  the  numerical  solu¬ 
tion  at  every  grid  point  in  the  spatial  domain  at  some  arbitrary  time  level. 
The  disturbance  is  expanded  into  a  Fourier  series  and  each  Fourier  component 
of  that  series  is  analyzed  separately.  The  FDE  is  stable  if  all  of  the 
Fourier  components  do  not  grow  in  time  and  unstable  if  any  one  of  the 
Fourier  components  grows  in  time.  The  method  of  analyzing  each  Fourier 
component  separately  is  valid  only  when  the  FDE's  are  linear  with  respect  to 
the  dependent  variable.  For  this  examination  the  FDE  is  not  linear  and  must 
be  linearized  before  the  stability  property  can  be  determined. 

A  Fourier  stability  analysis  has  been  performed  on  the  upwind  scheme 
using  the  linear  wave  equation.  Results  show  that  it  is  conditionally 
stable  with  the  following  conditions: 

0  <  j  =  c  At/A  x  <2  ( B.  8 ) 


4.  CONVERGENCE 

An  FDE  is  convergent  if  the  numerical  solution  of  the  FDE  approaches  the 
exact  solution  of  the  PDE  as  the  time-step  size  and  grid  spacing  approach 
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zero.  A  convergent  FDE  can  yield  a  solution  of  any  desired  accuracy  by 
reducing  the  time-step  size  (A  t)  and  the  grid  spacing  (A  x).  The  analysis 
of  the  convergence  of  an  FDE  for  a  complex  PDE  is  extremely  difficult,  and 
convergence  analysis  techniques  are  only  available  for  linear  PDE's. 

The  convergence  analysis  has  been  performed  for  the  upwind  scheme  with 
respect  to  the  linear  wave  equation  and  results  show  the  FDE  to  be  conver¬ 
gent  (Reference  1). 

If  the  FDE  is  well-posed,  consistent,  and  stable  (as  is  the  case  in  this 
examination);  the  Lax  Equivalence  Theorem  states  that  the  FDE  of  the  linear 
PDE  is  convergent.  Therefore,  according  to  this  theorem,  the  Warming-Beam 
upwind  scheme  should  be  convergent.  For  FDE's  of  quasi-linear  and  nonlinear 
PDE's,  the  Lax  Equivalence  Theorem  serves  as  an  important  guideline?  hence, 
consistency  and  stability  are  crucial  tests  for  convergence. 

5.  MODIFIED  EQUATION 

In  the  modified  equation  analysis,  a  PDE  is  derived  that  is  mathema¬ 
tically  equivalent  to  the  FDE  to  be  examined.  The  resulting  PDE  is  called 
the  modified  equation. 

The  modified  equation  is  derived  by  the  following  two-step  procedure: 

(1)  Expand  each  term  in  the  FDE  in  a  Taylor  series  expansion. 

(2)  Express  all  time  derivatives  (with  the  exception  of  the  first- 
order  time  derivative)  in  terms  of  spatial  derivatives. 

The  modified  equation  for  this  examination  is  given  by 

2  4  2 

U  +CU  =  (cAx  /6)  ( ] - v )  (2-v)  U  -  (  x  /8At) v  (1-v)  (2-v)  U  +...(B.9) 

t  X  XXX  xxxx 

This  method  of  analysis  is  used  exclusively  in  support  of  the  dissipation 
and  dispersion  error  investigations. 

6.  PHASE  AND  DISPERSION  ERROR 

Dispersion  is  mathematically  described  by  the  odd-order  spatial  deriva¬ 
tives.  The  coefficients  of  the  odd-order  spatial  derivatives  in  the 
modified  equation  for  the  FDE  must  be  identical  to  the  corresponding  coeffi¬ 
cients  in  the  PDE  in  order  for  a  FDE  to  have  the  same  dispersive  character¬ 
istics  as  the  PDE  it  is  to  represent.  If  the  corresponding  coefficients  are 
different,  then  the  solutions  of  the  FDE's  contain  dispersion  errors.  Since 
the  coefficients  of  the  odd-order  spatial  derivatives  in  the  modified  equa¬ 
tion  do  not  match  the  corresponding  coefficients  in  the  PDE,  the  FDE  con¬ 
tains  dispersion  error.  The  order  of  dispersion  is  equal  to  the  order  of 


the  coefficients  of  the  lowest  odd-order  spatial  derivative  excluding  the 
first  order  spatial  derivative.  Therefore,  the  FDE  contains  second-order 
dispersion;  the  higher  the  order  of  dispersion,  the  lower  the  dispersion 


U  +  CU  =0 
t  x 

U  +  CU  -  (cAx2/6)  ( 1-v)  (2-v)  U  +  (Ax4/ 8At)  vf1  v)2(2-v)  U  +•••  =  0 

t  X  XXX  xxxx 

The  dispersion  error  is  examined  by  the  Von  Neumann  method  (Reference 
1 ).  The  dispersion  error  for  each  Fourier  component  of  a  disturbance  after 
n  time  steps  is  given  by 


n(0pde  -  #f de ) 


(B. 10) 


where 


Opde  =  phase  angle  of  the  amplification  factor  ( PDE ) 
Ofde  =  phase  angle  of  the  amplification  factor  (FDE) 


The  dispersion  error  after  n  time  steps,  for  this  analysis,  is  given  by 


-'(  ■>  -  tan 


l-2v  (v+2  (1-v)  sin  2  )  sin  2 

2  Y 

vs in  y (1+2  (1-v)  sin  j ) 


(B.11 ) 


The  relative  phase  shift  error  for  a  given  Fourier  component  after  one 
time  step  is 


Ofde  /  Opde 


(B- 12) 


The  relative  phase  shift  error  for  this  examination  is 

-1  2  —  2  — 
tan  1-2  v(v+2  (1-v)  sin  2)  sin  2 

2  y 

vsin  y  (1+2  (1-v)  sin  -  ) 


(B. 13) 


The  dispersion  error  is  given  by  these  two  relationships  because  the  phase 
angle  of  the  amplification  factor  depends  only  on  the  odd-order  spatial 
derivatives  when  the  highest  time  derivative  is  first-order. 

For  leading  phase  error,  the  relative  phase  shift  error  must  be  greater 
than  unity  for  a  given  Fourier  component  (the  numerical  solution  for  that 
Fourier  component  gives  a  wave  speed  greater  than  the  wave  speed  given  by 
the  exact  solution).  Lagging  phase  error  results  when  the  relative  phase 
shift  error  is  less  than  unity  for  a  given  Fourier  component  (the  numerical 
solution  for  that  Fourier  component  gives  a  wave  speed  less  than  the  wave 
speed  given  by  the  exact  solution). 


The  Warming-Beam  upwind  scheme,  with  respect  to  the  convection  equation, 
has  a  lagging  phase  error  when  v  is  greater  than  one  and  a  predominantly 


leading  phase  error  when  is  less  than  one.  When  v  does  not  equal  one,  the 
relative  phase  shift  error  increases  as  the  wave  number,  y  (Kj  ax), 
increases  (Figure  B-1). 

7.  ARTIFICIAL  DISSIPATION 

As  diffusion  spreads  a  disturbance  in  every  direction,  the  disturbance 
is  smoothed  out  over  an  increasingly  large  area.  Diffusion  reduces  spatial 
gradients  and  lowers  the  magnitude  of  the  disturbance  by  spreading  it  out; 
this  phenomenon  is  termed  dissipation.  Diffusion  is  mathematically  des¬ 
cribed  by  even-order  spatial  derivatives  in  a  PDE;  if  the  even-order  spatial 
derivatives  in  the  PDE  are  zero,  the  PDE  will  not  have  any  dissipation. 

In  order  for  a  FDE  to  have  the  same  dissipative  characteristics  as  that 
of  the  PDE  it  is  modeling,  the  coefficients  of  the  even-order  spatial  deriv¬ 
atives  in  the  modified  equation  for  the  FDE  must  be  identical  to  the  corre¬ 
sponding  coefficients  of  the  PDE.  The  PDE  for  our  examination  is,  again, 

Ut  +  cux  =  0 

The  modified  equation  for  the  upwind  scheme  is  given  by  (Equationf  B.9) 

U  +  CU  =  (cAx2/6)  (1-v)  (2-v)  U  +  (Ax4/8At)v  (l-v)2(2-v)  U  +  •  •  • 
t  x  XXX  xxxx 

Therefore,  since  the  corresponding  even-order  coefficients  are  not  identical, 
there  is  fourth-order  dissipation  error  in  the  solution  of  the  FDE. 

The  dissipation  error  is  examined  using  the  Von  Neumann  method.  The 
dissipation  error  for  each  Fourier  component  of  a  disturbance  after  n  time 
steps  is 


(  Gr 


)  A°, 


Jpde 


Jfde 


pde  fde 

amplification  factor  (PDE) 

amplification  factor  (FDE) 

initial  amplitude  of  the  Fourier  component 


( B. 1 4 ) 


The  dissipation  error  is  given  by  this  relationship  because  the  modulus  of 
the  amplification  factor  depends  on  the  even-order  spatial  derivatives  when 
the  highest  time -derivative  is  first-order. 

The  dissipation  error,  for  this  examination,  for  a  given  time-step,  n, 
is  given  by 

1-  [  (1-4  v (1-  v)2  (2-v)  sin4  *)  1/2  ]n  (B.15) 

For  values  of  v  less  than  one,  as  decreases,  the  dissipation  error 
increases,  and  for  values  of  v  greater  than  one,  as  v  increases,  the  dis¬ 
sipation  error  inccreases.  Also,  for  values  of  v  not  equal  to  unity,  as 
the  wave  number,  y  (k.  Ax),  increases  the  dissipation  error  increases 
(Figure  B-2).  3 
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Figure  B-l.  Relative  Phase  Shift  Error 
vs  Courant  Number  ( m) 
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Figure 


AMPLIFICATION  FACTOR: 


B-2.  Amplification  Factor  Modulus 
vs  Courant  Number  (v) 
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